## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6     ✔ purrr   0.3.4
## ✔ tibble  3.2.1     ✔ dplyr   1.1.2
## ✔ tidyr   1.2.0     ✔ stringr 1.4.1
## ✔ readr   2.1.2     ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## 
## Attaching package: 'lubridate'
## 
## 
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
## 
## 
## corrplot 0.92 loaded
## 
## 
## Attaching package: 'ggpp'
## 
## 
## The following object is masked from 'package:ggplot2':
## 
##     annotate
## 
## 
## 
## Attaching package: 'flextable'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     compose
## 
## 
## 
## Attaching package: 'maps'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     map
## 
## 
## ℹ Google's Terms of Service: ]8;;https://mapsplatform.google.com<https://mapsplatform.google.com>]8;;
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
## 
## Attaching package: 'scales'
## 
## 
## The following object is masked from 'package:purrr':
## 
##     discard
## 
## 
## The following object is masked from 'package:readr':
## 
##     col_factor
## 
## 
## Loading required package: sp
## 
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired by the end of 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
## 
## 
## Attaching package: 'raster'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## rgeos version: 0.6-3, (SVN revision 696)
##  GEOS runtime version: 3.10.2-CAPI-1.16.0 
##  Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
##  GEOS using OverlayNG
##  Linking to sp version: 1.6-0 
##  Polygon checking: TRUE 
## 
## 
## 
## Attaching package: 'rgeos'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     symdiff

Overview

This script reads in cleaned datafiles with stream temperature and invertebrate data produced by the invert_data_cleaning.Rmd and teton_data_cleaning.Rmd scripts and performs more involved analyses of stream temperature and invertebrate density/diversity over time. A few findings of interest from previous analyses are that:

To build on these findings, the main sections of this script are:

I. Temperature:

  1. Trends in summer mean, max, min, and degree-days over time: Site-by-site and grouped by source

  2. Modeling stream temperature using SNOTEL data: Site-by-site and grouped by source

  3. Linear mixed models of temperature over time

II. Invertebrates:

  1. Trends in richness, abundance, and density over time: Site-by-site and grouped by source

  2. Modeling richness, abundance, and density using SNOTEL data: Site-by-site and grouped by source

III. Temperature and Invertebrates

  1. Test for correlations between temperature metrics (mean/max/min summer; degree days): Site-by-site and grouped by source

Read in necessary files:

stream_temp <- read.csv("Temperature/cleaned_full_datasets/temps_hourly.csv") #hourly stream temp data
invert_density <- read.csv("invert_data/cleaned_csv/full_invert_densities.csv")%>%
  rename(site = Stream)#invert density data
invert_richness <- read.csv("invert_data/cleaned_csv/invert_richness.csv") #invert richness data
snotel <- read.csv("teton_snotel.csv") #snotel data
sources <- read.csv("source_info.csv")%>%
  rename(site = stream) #source info, rename for merge

### adding source info to stream temps and invert datasets: 

stream_source <- merge(stream_temp, sources, all = T)%>%
  mutate(date1 = ymd(date1), 
         year = year(date1))%>%
  filter(!is.na(temp_c))#Add source info to hourly temperature data; re-code date as r date format
stream_source_daily <- stream_source%>%
  group_by(site, full_name, stream_code, date1, year, lat, long)%>%
  summarise(t_xbar = mean(temp_c, na.rm = T), 
            t_max = max(temp_c, na.rm = T), 
            t_min = min(temp_c, na.rm = T), 
            t_range = t_max-t_min, 
            source = unique(source)) #calculate daily mean, min, and max temps for each site along with temp range
## `summarise()` has grouped output by 'site', 'full_name', 'stream_code',
## 'date1', 'year', 'lat'. You can override using the `.groups` argument.
density_source <- merge(invert_density, sources, all = T)
richness_source <- merge(invert_richness, sources, all = T)

I. Temperature:

2. Modeling stream temperature using SNOTEL data

First, calculate April 1 SWE and summer air temp from snotel data; add to mean temperature dataset:

#Calc means across both sites:
snotel_means <- snotel %>%
  mutate(date1 = as.Date(date, format = "%m/%d/%y"), #recode date column as R date
         airtemp_c = as.numeric(airtemp_c))%>% #recode airtemp as numeric
  group_by(date1)%>% #group by date
  summarise(swe_xbar = mean(swe_cm), #calculate mean SWE and airtemp across both stations
            airtemp_xbar = mean(airtemp_c, na.rm = T))%>%
  mutate(mo = month(date1), yr = year(date1)) #Add month and year cols
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `airtemp_c = as.numeric(airtemp_c)`.
## Caused by warning:
## ! NAs introduced by coercion
#Calculate max SWE and mean summer airtemp for each year:
snotel_summary <- snotel_means %>% 
  group_by(yr)%>% #group by year
  mutate(swe_max = max(swe_xbar), 
         swe_max_date = date1[which.max(swe_xbar)], 
         swe_zero_date = date1[which.min(swe_xbar)][1])%>% #extract max SWE and add to new col
  ungroup()%>% #ungroup to get all columns back
  filter(mo == 7 | mo == 8)%>% #extract summer months
  group_by(yr)%>% #group by year
  summarise(airtemp_summer = mean(airtemp_xbar), #calculate mean summer air temp
            swe_max = unique(swe_max), 
            swe_max_date = ymd(unique(swe_max_date)), 
            swe_zero_date =ymd(unique(swe_zero_date)))%>% #keep max SWE
  mutate(melt_days = as.numeric(swe_zero_date-swe_max_date))%>%
  rename(year = yr)

#Merge with long-term temperature averages: 
stream_snotel <- merge(summer_3plus, snotel_summary)

A. Site-by-site

i. Correlations between SWE and Summer Temp:

LMs - first, nest by site:

stream_snotel_nested <- stream_snotel %>%
  nest(data = - site)

lm for each site:

swe.temp.lms <- stream_snotel_nested %>%
  mutate(model = purrr::map(data, ~lm(summer_xbar~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "swe_max")  #remove intercept term to just get info on year for each site
swe.temp.lms
## # A tibble: 12 × 7
##    site       data              term     estimate std.error statistic p.value
##    <chr>      <list>            <chr>       <dbl>     <dbl>     <dbl>   <dbl>
##  1 s_cascade  <tibble [5 × 18]> swe_max -0.0185     0.0454     -0.407  0.711 
##  2 n_teton    <tibble [8 × 18]> swe_max -0.0455     0.0332     -1.37   0.219 
##  3 mid_teton  <tibble [8 × 18]> swe_max  0.000589   0.00313     0.188  0.857 
##  4 windcave   <tibble [5 × 18]> swe_max -0.00725    0.00771    -0.941  0.416 
##  5 s_ak_basin <tibble [6 × 18]> swe_max -0.00762    0.0101     -0.753  0.494 
##  6 paintbrush <tibble [6 × 18]> swe_max -0.0477     0.0400     -1.19   0.299 
##  7 delta      <tibble [6 × 18]> swe_max -0.00137    0.00313    -0.438  0.684 
##  8 silver_run <tibble [5 × 18]> swe_max -0.0324     0.00963    -3.37   0.0435
##  9 s_teton    <tibble [6 × 18]> swe_max -0.0186     0.0311     -0.598  0.582 
## 10 grizzly    <tibble [5 × 18]> swe_max -0.0395     0.0488     -0.810  0.477 
## 11 skillet    <tibble [5 × 18]> swe_max -0.0211     0.0174     -1.22   0.311 
## 12 cloudveil  <tibble [4 × 18]> swe_max -0.0174     0.0191     -0.913  0.457

Add P values and merge:

swe_summer_ps <- merge(swe.temp.lms, stream_snotel)

swe.summer.temp <- ggplot(swe_summer_ps, aes(x = swe_max, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Annual maxium SWE, cm", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'

ii. Summer stream temp ~ summer air temp:

lm for each site:

air.temp.lms <- stream_snotel_nested %>%
  mutate(model = purrr::map(data, ~lm(summer_xbar~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "airtemp_summer")  #remove intercept term to just get info on year for each site
air.temp.lms
## # A tibble: 12 × 7
##    site       data              term        estimate std.error statistic p.value
##    <chr>      <list>            <chr>          <dbl>     <dbl>     <dbl>   <dbl>
##  1 s_cascade  <tibble [5 × 18]> airtemp_su… -0.683      0.270    -2.53    0.0852
##  2 n_teton    <tibble [8 × 18]> airtemp_su… -0.880      0.815    -1.08    0.322 
##  3 mid_teton  <tibble [8 × 18]> airtemp_su… -0.0492     0.0708   -0.695   0.513 
##  4 windcave   <tibble [5 × 18]> airtemp_su…  0.00876    0.240     0.0365  0.973 
##  5 s_ak_basin <tibble [6 × 18]> airtemp_su… -0.405      0.305    -1.33    0.255 
##  6 paintbrush <tibble [6 × 18]> airtemp_su…  0.110      1.73      0.0637  0.952 
##  7 delta      <tibble [6 × 18]> airtemp_su…  0.112      0.105     1.07    0.346 
##  8 silver_run <tibble [5 × 18]> airtemp_su… -0.929      1.59     -0.585   0.600 
##  9 s_teton    <tibble [6 × 18]> airtemp_su… -0.123      1.20     -0.103   0.923 
## 10 grizzly    <tibble [5 × 18]> airtemp_su… -0.339      1.42     -0.238   0.827 
## 11 skillet    <tibble [5 × 18]> airtemp_su…  0.584      0.456     1.28    0.290 
## 12 cloudveil  <tibble [4 × 18]> airtemp_su… -0.224      0.451    -0.497   0.668

Add P values and plot:

air_summer_ps <- merge(air.temp.lms, stream_snotel)

swe.summer.temp <- ggplot(air_summer_ps, aes(x = airtemp_summer, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
swe.summer.temp
## `geom_smooth()` using formula 'y ~ x'

iii. Summer stream temp ~ meltout length (days):

lm for each site:

melt.temp.lms <- stream_snotel_nested %>%
  mutate(model = purrr::map(data, ~lm(summer_xbar~melt_days, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "melt_days")  #remove intercept term to just get info on year for each site
melt.temp.lms
## # A tibble: 12 × 7
##    site       data              term       estimate std.error statistic p.value
##    <chr>      <list>            <chr>         <dbl>     <dbl>     <dbl>   <dbl>
##  1 s_cascade  <tibble [5 × 18]> melt_days  0.00638    0.0768     0.0831  0.939 
##  2 n_teton    <tibble [8 × 18]> melt_days -0.144      0.0537    -2.68    0.0364
##  3 mid_teton  <tibble [8 × 18]> melt_days -0.00883    0.00550   -1.61    0.160 
##  4 windcave   <tibble [5 × 18]> melt_days -0.0379     0.0185    -2.04    0.134 
##  5 s_ak_basin <tibble [6 × 18]> melt_days  0.0308     0.0127     2.43    0.0721
##  6 paintbrush <tibble [6 × 18]> melt_days -0.0532     0.0828    -0.642   0.556 
##  7 delta      <tibble [6 × 18]> melt_days  0.000493   0.00599    0.0824  0.938 
##  8 silver_run <tibble [5 × 18]> melt_days -0.0251     0.0296    -0.849   0.458 
##  9 s_teton    <tibble [6 × 18]> melt_days -0.101      0.0332    -3.05    0.0381
## 10 grizzly    <tibble [5 × 18]> melt_days -0.0988     0.0476    -2.08    0.129 
## 11 skillet    <tibble [5 × 18]> melt_days -0.0184     0.0273    -0.675   0.548 
## 12 cloudveil  <tibble [4 × 18]> melt_days -0.00711    0.0242    -0.293   0.797

Add P values and plot:

melt_summer_ps <- merge(melt.temp.lms, stream_snotel)

melt.summer.temp <- ggplot(melt_summer_ps, aes(x = melt_days, y = summer_xbar, color = source))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~full_name, scales = "free")+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  geom_text_npc(aes(npcx = 0.9, npcy = 0.9, label = paste( "P = ",round(p.value, 2), sep = "")))+
  labs(x = "Mean daily July-Aug. air temp, C", y = "Mean daily July-Aug. stream temp., C", color = "Stream Source")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))
melt.summer.temp
## `geom_smooth()` using formula 'y ~ x'

B. Grouped by source

First, add SNOTEL summary info to mean summer temps grouped by source:

source_snotel <- merge(summer_source, snotel_summary)
i. LMs of mean/min/max temp ~ SWE:

Re-organize data for analysis and plotting:

source_swe <- source_snotel %>%
  dplyr::select(source, swe_max, t_xbar, t_min, t_max, t_range)%>% #dplyr::select temperature, source, and swe cols
  pivot_longer(!c(source, swe_max), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

Nest data by source and measure:

sourceSwe_nested <- source_swe %>%
  nest(data = -c(source, measure))

LMs for each source and measure:

source.swe.lms <- sourceSwe_nested %>%
  mutate(model = purrr::map(data, ~lm(temp~swe_max, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "swe_max")  #remove intercept term to just get info on year for each site
source.swe.lms
## # A tibble: 12 × 8
##    source           measure data     term   estimate std.error statistic p.value
##    <chr>            <chr>   <list>   <chr>     <dbl>     <dbl>     <dbl>   <dbl>
##  1 Glacier          t_xbar  <tibble> swe_m… -0.00443   0.00821    -0.540  0.609 
##  2 Glacier          t_min   <tibble> swe_m… -0.00477   0.00763    -0.626  0.554 
##  3 Glacier          t_max   <tibble> swe_m… -0.00108   0.00992    -0.109  0.917 
##  4 Glacier          t_range <tibble> swe_m…  0.00369   0.00767     0.481  0.647 
##  5 Snow melt        t_xbar  <tibble> swe_m… -0.0430    0.0250     -1.72   0.136 
##  6 Snow melt        t_min   <tibble> swe_m… -0.0434    0.0191     -2.28   0.0629
##  7 Snow melt        t_max   <tibble> swe_m… -0.0354    0.0331     -1.07   0.326 
##  8 Snow melt        t_range <tibble> swe_m…  0.00807   0.0188      0.429  0.683 
##  9 Subterranean ice t_xbar  <tibble> swe_m…  0.00273   0.0131      0.208  0.844 
## 10 Subterranean ice t_min   <tibble> swe_m…  0.00956   0.00674     1.42   0.215 
## 11 Subterranean ice t_max   <tibble> swe_m… -0.0101    0.0276     -0.366  0.729 
## 12 Subterranean ice t_range <tibble> swe_m… -0.0197    0.0242     -0.812  0.454

Add p-values and plot:

sourceSwe_pval<-merge(source_swe, source.swe.lms)%>%
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85, measure == "t_range"~0.8)) #create column for adjusting P-value label Y position

source.swe<- ggplot(sourceSwe_pval, aes(x = swe_max, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkorange","darkgreen"), labels = c("Max. temp.", "Min. temp.", "Temp. range", "Mean temp"))+
  labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.9, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.swe
## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_swe.pdf", source.swe, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
ii. LMs of mean/min/max temp ~ Air temp:

Re-organize data for analysis and plotting:

source_air <- source_snotel %>%
  dplyr::select(source, airtemp_summer, t_xbar, t_min, t_max, t_range)%>% #dplyr::select temperature, source, and swe cols
  pivot_longer(!c(source, airtemp_summer), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

Nest data by source and measure:

sourceAir_nested <- source_air %>%
  nest(data = -c(source, measure))

LMs for each source and measure:

source.air.lms <- sourceAir_nested %>%
  mutate(model = purrr::map(data, ~lm(temp~airtemp_summer, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "airtemp_summer")  #remove intercept term to just get info on year for each site
source.air.lms
## # A tibble: 12 × 8
##    source           measure data     term   estimate std.error statistic p.value
##    <chr>            <chr>   <list>   <chr>     <dbl>     <dbl>     <dbl>   <dbl>
##  1 Glacier          t_xbar  <tibble> airte…    0.208     0.178     1.17   0.287 
##  2 Glacier          t_min   <tibble> airte…    0.240     0.156     1.53   0.176 
##  3 Glacier          t_max   <tibble> airte…    0.103     0.229     0.451  0.668 
##  4 Glacier          t_range <tibble> airte…   -0.137     0.174    -0.784  0.463 
##  5 Snow melt        t_xbar  <tibble> airte…   -1.11      0.554    -2.01   0.0917
##  6 Snow melt        t_min   <tibble> airte…   -0.955     0.470    -2.03   0.0882
##  7 Snow melt        t_max   <tibble> airte…   -1.19      0.694    -1.71   0.138 
##  8 Snow melt        t_range <tibble> airte…   -0.232     0.437    -0.532  0.614 
##  9 Subterranean ice t_xbar  <tibble> airte…   -0.525     0.169    -3.10   0.0270
## 10 Subterranean ice t_min   <tibble> airte…   -0.156     0.161    -0.970  0.377 
## 11 Subterranean ice t_max   <tibble> airte…   -1.15      0.337    -3.42   0.0189
## 12 Subterranean ice t_range <tibble> airte…   -0.995     0.350    -2.84   0.0363

Add p-values and plot:

sourceAir_pval<-merge(source_air, source.air.lms)%>%
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85, measure == "t_range"~0.8)) #create column for adjusting P-value label Y position

source.air<- ggplot(sourceAir_pval, aes(x = airtemp_summer, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkorange","darkgreen"), labels = c("Max. temp.", "Min. temp.", "Temp. range", "Mean temp"))+
  labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
  labs(x = "Mean July-Aug. air temp., C", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.air
## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_air.pdf", source.air, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'
iii. LMs of mean/min/max temp ~ Meltout length:

Re-organize data for analysis and plotting:

source_melt <- source_snotel %>%
  dplyr::select(source, melt_days, t_xbar, t_min, t_max, t_range)%>% #dplyr::select temperature, source, and swe cols
  pivot_longer(!c(source, melt_days), names_to = "measure", values_to = "temp") #pivot longer for plotting + analysis

Nest data by source and measure:

sourceMelt_nested <- source_melt %>%
  nest(data = -c(source, measure))

LMs for each source and measure:

source.melt.lms <- sourceMelt_nested %>%
  mutate(model = purrr::map(data, ~lm(temp~melt_days, data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     tidy))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model) %>% #unnest model column to get rows for each site, columns for model metrics
  filter(term == "melt_days")  #remove intercept term to just get info on year for each site
source.melt.lms
## # A tibble: 12 × 8
##    source           measure data     term   estimate std.error statistic p.value
##    <chr>            <chr>   <list>   <chr>     <dbl>     <dbl>     <dbl>   <dbl>
##  1 Glacier          t_xbar  <tibble> melt_… -0.00818    0.0173   -0.473    0.653
##  2 Glacier          t_min   <tibble> melt_… -0.0123     0.0157   -0.786    0.462
##  3 Glacier          t_max   <tibble> melt_… -0.00128    0.0208   -0.0614   0.953
##  4 Glacier          t_range <tibble> melt_…  0.0111     0.0157    0.703    0.508
##  5 Snow melt        t_xbar  <tibble> melt_… -0.0751     0.0561   -1.34     0.230
##  6 Snow melt        t_min   <tibble> melt_… -0.0560     0.0495   -1.13     0.301
##  7 Snow melt        t_max   <tibble> melt_… -0.0937     0.0652   -1.44     0.201
##  8 Snow melt        t_range <tibble> melt_… -0.0377     0.0369   -1.02     0.347
##  9 Subterranean ice t_xbar  <tibble> melt_…  0.0138     0.0252    0.546    0.609
## 10 Subterranean ice t_min   <tibble> melt_…  0.00267    0.0157    0.170    0.871
## 11 Subterranean ice t_max   <tibble> melt_…  0.0290     0.0536    0.540    0.612
## 12 Subterranean ice t_range <tibble> melt_…  0.0263     0.0494    0.532    0.618

Add p-values and plot:

sourceMelt_pval<-merge(source_melt, source.melt.lms)%>%
  mutate(ypos = case_when(measure == "t_xbar" ~ 0.9, measure == "t_max" ~ 0.95, measure == "t_min" ~ 0.85, measure == "t_range"~0.8)) #create column for adjusting P-value label Y position

source.melt<- ggplot(sourceMelt_pval, aes(x = melt_days, y = temp, color = measure))+
  geom_point()+
  geom_line()+
  geom_smooth(se = F, method = "lm", linetype = "dashed", size = 0.5)+
  theme_classic()+
  facet_wrap(~source, scales = "free_x")+
   scale_color_manual(values = c("darkred","darkblue", "darkorange","darkgreen"), labels = c("Max. temp.", "Min. temp.", "Temp. range", "Mean temp"))+
  labs(x = "Annual maximum SWE, cm", y = "Stream temperature, C", color = "Measurement Type")+
  labs(x = "Days from peak SWE to zero SWE", y = "Stream temperature, C", color = "Measurement Type")+
  theme(legend.position = "bottom", 
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1), 
        axis.title = element_text(size = 11, face = "bold"), 
        axis.text = element_text(size = 8), 
        legend.text = element_text(size = 10), 
        legend.title = element_text(size = 11, face = "bold"), 
        strip.text = element_text(size = 10, face = "bold"))+
  geom_text_npc(aes(npcx = 0.1, npcy = ypos, color = measure, label = paste( "P = ",round(p.value, 3), sep = ""))) # add p value labels to plot
source.melt
## `geom_smooth()` using formula 'y ~ x'

Save:

ggsave("Temperature/plots/source_meltout.pdf", source.melt, device = "pdf", units = "in", height = 4, width = 6.5, dpi = "retina")
## `geom_smooth()` using formula 'y ~ x'

3. Linear mixed models of temperature over time

General concept: want to test a linear mixed model of stream temperature (min, max, mean, and range) with:

  • Fixed effects:
    • Year
    • Max SWE
    • Max SWE day of year
    • Mean summer air temp
    • Meltout length (days)
  • Random intercept for Site

Current approach will be to run seperate models for each different source type.

Calculate day of year for max SWE and zero SWE:

stream_snotel1 <- stream_snotel%>%
  mutate(swe_max_doy =yday(swe_max_date), 
         swe_zero_doy=yday(swe_zero_date))

A. Testing for correlation between predictor variables:

Extract predictor variables:

predictors <- stream_snotel1%>%
  dplyr::select(year, airtemp_summer, swe_max, melt_days, swe_max_doy, swe_zero_doy)%>%
  distinct()

Calculate correlation matrix and melt for plotting:

cormat <-round(cor(predictors), 2)
cormat[upper.tri(cormat)]<- NA
cormat_melt<-reshape2::melt(cormat, na.rm = T)
cormat_melt
##              Var1           Var2 value
## 1            year           year  1.00
## 2  airtemp_summer           year  0.75
## 3         swe_max           year -0.18
## 4       melt_days           year -0.32
## 5     swe_max_doy           year  0.58
## 6    swe_zero_doy           year  0.22
## 8  airtemp_summer airtemp_summer  1.00
## 9         swe_max airtemp_summer  0.20
## 10      melt_days airtemp_summer -0.28
## 11    swe_max_doy airtemp_summer  0.72
## 12   swe_zero_doy airtemp_summer  0.42
## 15        swe_max        swe_max  1.00
## 16      melt_days        swe_max  0.20
## 17    swe_max_doy        swe_max  0.26
## 18   swe_zero_doy        swe_max  0.55
## 22      melt_days      melt_days  1.00
## 23    swe_max_doy      melt_days -0.68
## 24   swe_zero_doy      melt_days  0.59
## 29    swe_max_doy    swe_max_doy  1.00
## 30   swe_zero_doy    swe_max_doy  0.19
## 36   swe_zero_doy   swe_zero_doy  1.00

Plot correlation matrix:

ggplot(cormat_melt, aes(x = Var1, y = Var2, fill = value))+
  geom_tile(color = "white")+
 scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
   midpoint = 0, limit = c(-1,1), space = "Lab", 
   name="Pearson\nCorrelation") +
  theme_minimal()+ 
 theme(axis.text.x = element_text(angle = 45, vjust = 1, 
    size = 10, hjust = 1), 
    axis.text.y = element_text(
    size = 10),
    axis.title = element_blank())+
 coord_fixed()+
  geom_text(aes(label = value))

Using a threshold of 0.5, correlations are:

  • melt_days~swe_max_doy
  • melt_days~swe_zero_doy
  • swe_zero_doy~swe_max
  • airtemp_summer~swe_max_doy
  • airtemp_summer~year
  • swe_max_doy~year

Problematic variables are swe max/zero DOY, year. For now, dropping year (more relevant info contained in summer air temperature) and Swe max/zero DOY (info from both is roughly contained in melt_days)

Confirming that remaining variables don’t have correlations greater than 0.5:

cormat_filter <- cormat_melt%>%
  filter(!(Var1 == "year"|Var1=="swe_max_doy"|Var1=="swe_zero_doy"),
         !(Var2 == "year"|Var2=="swe_max_doy"|Var2=="swe_zero_doy"))

ggplot(cormat_filter, aes(x = Var1, y = Var2, fill = value))+
  geom_tile(color = "white")+
 scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
   midpoint = 0, limit = c(-1,1), space = "Lab", 
   name="Pearson\nCorrelation") +
  theme_minimal()+ 
 theme(axis.text.x = element_text(angle = 45, vjust = 1, 
    size = 10, hjust = 1), 
    axis.text.y = element_text(
    size = 10),
    axis.title = element_blank())+
 coord_fixed()+
  geom_text(aes(label = value))

B. Checking distribution of the response variables

Data prep for plotting & modeling:

lmm_pivot <- stream_snotel1 %>%
  dplyr::select(max_xbar, min_xbar, range_xbar, summer_xbar, source, airtemp_summer, swe_max, melt_days, site, year)%>%
  pivot_longer(!c(source, airtemp_summer, swe_max, melt_days, site, year), names_to = "measure", values_to = "temp_c")%>%
  mutate(transform = log10(temp_c))

#for plotting model results: 

lmm_plot <- lmm_pivot%>%
  pivot_longer(!c(source, site, year, measure, temp_c, transform), names_to = "predictor_name", values_to = "predictor_value")

Check distribution of response variables:

#without log10 transform:
ggplot(lmm_pivot, aes(x = temp_c))+
  geom_density(fill = "grey", color = "black")+
  facet_wrap(~measure)+
  theme_classic()

#with log10 transform:
ggplot(lmm_pivot, aes(x = transform))+
  geom_density(fill = "grey", color = "black")+
  facet_wrap(~measure)+
  theme_classic()

All are left-skewed, but a log transform just makes them left-skewed. Skipping log transform for now.

C. Running LMMs

Load packages:

library(lme4)
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## Attaching package: 'lme4'
## The following object is masked from 'package:raster':
## 
##     getData
library(lmerTest)
## 
## Attaching package: 'lmerTest'
## The following object is masked from 'package:lme4':
## 
##     lmer
## The following object is masked from 'package:stats':
## 
##     step

Nest data by source:

lmm_nested <- lmm_pivot%>%
  nest(data = -c(source, measure))
i. Maximum temp:

Run LMM for each source and measure - form:

(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)

max_temps <- filter(lmm_nested, measure == "max_xbar")

max.lmm <- max_temps %>%
  mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
  mutate(p.value = round(p.value, 4))%>%
  filter(!term == "(Intercept)")%>%
  dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
## boundary (singular) fit: see help('isSingular')
max.lmm
## # A tibble: 9 × 7
##   source   term           estimate std.error statistic    df p.value
##   <chr>    <chr>             <dbl>     <dbl>     <dbl> <dbl>   <dbl>
## 1 sub_ice  melt_days      -0.0138    0.0421    -0.328   12.0  0.748 
## 2 sub_ice  airtemp_summer -1.18      0.429     -2.74    12.0  0.0179
## 3 sub_ice  swe_max         0.00178   0.0184     0.0968  12.0  0.924 
## 4 snowmelt melt_days      -0.120     0.0375    -3.19    22.2  0.0042
## 5 snowmelt airtemp_summer -1.70      0.620     -2.74    22.6  0.0117
## 6 snowmelt swe_max        -0.0219    0.0199    -1.10    22.5  0.282 
## 7 glacier  melt_days      -0.00648   0.0155    -0.417   16.0  0.682 
## 8 glacier  airtemp_summer -0.0621    0.226     -0.275   16.1  0.787 
## 9 glacier  swe_max         0.00214   0.00886    0.242   16.1  0.812

Results:

Subterranean ice:

Significant negative correlation between summer air temp and max stream temp (maybe more ice melting?)

Snowmelt:

Significant negative correlation between melt length and max temp; same with summer air temp and max temp (??)

Glacier:

No significant effects.

Visualization:

ggplot(filter(lmm_plot, measure == "max_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
  geom_point(size = 0.5)+
  geom_line(size = 0.5)+
  geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
  theme_classic()+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  facet_grid(source~predictor_name, scales = "free")+
  theme(legend.position = "none")+
  labs(x = "Predictor Value", y = "Mean maximum July-Aug stream temp, C")
## `geom_smooth()` using formula 'y ~ x'

ii. Minimum temp:

Run LMM for each source and measure - form:

(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)

min_temps <- filter(lmm_nested, measure == "min_xbar")

min.lmm <- min_temps %>%
  mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
  mutate(p.value = round(p.value, 4))%>%
  filter(!term == "(Intercept)")%>%
  dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
min.lmm
## # A tibble: 9 × 7
##   source   term            estimate std.error statistic    df p.value
##   <chr>    <chr>              <dbl>     <dbl>     <dbl> <dbl>   <dbl>
## 1 sub_ice  melt_days      -0.000123   0.00849   -0.0144  10.0  0.989 
## 2 sub_ice  airtemp_summer -0.0525     0.0888    -0.591   10.0  0.568 
## 3 sub_ice  swe_max        -0.00722    0.00376   -1.92    10.0  0.0838
## 4 snowmelt melt_days      -0.0684     0.0128    -5.34    22.2  0     
## 5 snowmelt airtemp_summer -0.750      0.212     -3.54    22.5  0.0018
## 6 snowmelt swe_max        -0.0221     0.00680   -3.24    22.5  0.0037
## 7 glacier  melt_days      -0.00525    0.00642   -0.818   16.0  0.426 
## 8 glacier  airtemp_summer  0.0540     0.0933     0.579   16.1  0.571 
## 9 glacier  swe_max        -0.00654    0.00366   -1.79    16.1  0.0927

Results:

Subterranean ice:

Significant negative correlation between maximum SWE and min summer stream temp (Maybe less ice melting bc of snow insulation?)

Snowmelt:

Significant negative correlation between melt length, summer air temp, and swe max and summer stream temp. Melt length makes sense, others don’t really…

Glacier:

No significant effects.

Visualization:

ggplot(filter(lmm_plot, measure == "min_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
  geom_point(size = 0.5)+
  geom_line(size = 0.5)+
  geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
  theme_classic()+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  facet_grid(source~predictor_name, scales = "free")+
  theme(legend.position = "none")+
  labs(x = "Predictor Value", y = "Mean minimum July-Aug stream temp, C")
## `geom_smooth()` using formula 'y ~ x'

iii. Average temp:

Run LMM for each source and measure - form:

(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)

mean_temps <- filter(lmm_nested, measure == "summer_xbar")

mean.lmm <- mean_temps %>%
  mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
  mutate(p.value = round(p.value, 4))%>%
  filter(!term == "(Intercept)")%>%
  dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
mean.lmm
## # A tibble: 9 × 7
##   source   term            estimate std.error statistic    df p.value
##   <chr>    <chr>              <dbl>     <dbl>     <dbl> <dbl>   <dbl>
## 1 sub_ice  melt_days      -0.000173   0.0187   -0.00927  10.1  0.993 
## 2 sub_ice  airtemp_summer -0.416      0.195    -2.13     10.1  0.0586
## 3 sub_ice  swe_max        -0.00406    0.00828  -0.491    10.1  0.634 
## 4 snowmelt melt_days      -0.0950     0.0215   -4.43     22.1  0.0002
## 5 snowmelt airtemp_summer -1.19       0.349    -3.41     24.2  0.0023
## 6 snowmelt swe_max        -0.0219     0.0112   -1.95     23.5  0.0635
## 7 glacier  melt_days      -0.00520    0.00863  -0.602    16.0  0.556 
## 8 glacier  airtemp_summer  0.0124     0.126     0.0989   16.1  0.922 
## 9 glacier  swe_max        -0.00429    0.00492  -0.871    16.1  0.397

Results:

Subterranean ice:

Significant negative correlation between summer air temp and mean summer stream temp (Maybe less ice melting bc of snow insulation?)

Snowmelt:

Significant negative correlation between melt length, summer air temp, and swe max and summer stream temp. Melt length makes sense, others don’t really…

Glacier:

No significant effects.

Visualization:

ggplot(filter(lmm_plot, measure == "summer_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
  geom_point(size = 0.5)+
  geom_line(size = 0.5)+
  geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
  theme_classic()+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  facet_grid(source~predictor_name, scales = "free")+
  theme(legend.position = "none")+
  labs(x = "Predictor Value", y = "Mean daily July-Aug stream temp, C")
## `geom_smooth()` using formula 'y ~ x'

iv. temp range:

Run LMM for each source and measure - form:

(t_max, t_min, t_range, t_ave) ~ airtemp_summer + swe_max + melt_days + random(site)

temp_ranges <- filter(lmm_nested, measure == "range_xbar")

range.lmm <- temp_ranges %>%
  mutate(model = purrr::map(data, ~lmerTest::lmer(temp_c~melt_days+airtemp_summer+swe_max+(1|site), data = .)%>% #run a model of temp~year for each site (referenced by "."); store results in new column called "model"
                     broom.mixed::tidy(effects="fixed", conf.int = TRUE)))%>% #tidy function to construct a summary table with model results (seperate tables for each site)
  unnest(model)%>% #unnest model column to get rows for each site, columns for model metrics
  mutate(p.value = round(p.value, 4))%>%
  filter(!term == "(Intercept)")%>%
  dplyr::select(source, term, estimate, std.error, statistic, df, p.value)
range.lmm
## # A tibble: 9 × 7
##   source   term           estimate std.error statistic    df p.value
##   <chr>    <chr>             <dbl>     <dbl>     <dbl> <dbl>   <dbl>
## 1 sub_ice  melt_days      -0.00424   0.0432    -0.0982  10.2  0.924 
## 2 sub_ice  airtemp_summer -0.977     0.450     -2.17    10.5  0.0538
## 3 sub_ice  swe_max         0.00448   0.0191     0.235   10.4  0.819 
## 4 snowmelt melt_days      -0.0508    0.0270    -1.89    22.1  0.0726
## 5 snowmelt airtemp_summer -0.977     0.447     -2.19    22.2  0.0397
## 6 snowmelt swe_max        -0.00179   0.0143    -0.125   22.2  0.902 
## 7 glacier  melt_days      -0.00113   0.0130    -0.0874  16.0  0.931 
## 8 glacier  airtemp_summer -0.114     0.189     -0.604   16.0  0.554 
## 9 glacier  swe_max         0.00861   0.00739    1.16    16.0  0.261

Results:

Subterranean ice:

Marginally significant negative correlation between summer air temp and temperature variability (warmer summers = less variability)

Snowmelt:

Significant negative correlation of melt length and summer air temp with temperature variability (longer melt and warmer summers = less variability). Seems a bit contradictory.

Glacier:

No significant effects.

Visualization:

ggplot(filter(lmm_plot, measure == "range_xbar"), aes(x = predictor_value, y = temp_c, color = source, group = site))+
  geom_point(size = 0.5)+
  geom_line(size = 0.5)+
  geom_smooth(aes(group = NULL), method = "lm", se = T, linetype = "dashed", color = "black", size = 0.5)+
  theme_classic()+
  scale_color_manual(values = pres.pal, labels = c("Glacier", "Snowmelt", "Icy Seep"))+
  facet_grid(source~predictor_name, scales = "free")+
  theme(legend.position = "none")+
  labs(x = "Predictor Value", y = "Mean July-Aug daily stream temp. change, C")
## `geom_smooth()` using formula 'y ~ x'

II. Invertebrates:

2. Modeling richness, abundance, and density using SNOTEL data

A. Site-by-site

B. Grouped by source

III. Temperature and Invertebrates

1. Test for correlations between temperature metrics (mean/max/min summer; degree days)

A. Site-by-site

B. Grouped by source